### Replication files for
### "When Democracy Brings Insecurity: The Political Legacies of Regime Change"
### Hans Lueders
### April 2023
### Part 1: Survey experiment



# Preparation -------------------------------------------------------------

### packages
library(ggplot2)
library(stargazer)
library(lfe)
library(xtable)
library(readxl)
library(tm)
library(stringr)
library(SnowballC)
library(rgdal)
library(MASS)



### set working directory
path <- "" # define path here
setwd(path)

### empty environment
rm(list=ls())


### load data
load("Data/WhenDemocracyBringsInsecurity_HansLueders_SurveyExperiment.RData")




# Figure 2 ----------------------------------------------------------------

### prepare output container
out <- data.frame(matrix(nrow = 9, ncol=4))
colnames(out) <- c("region", "treatment", "estimate", "se")
out$region <- rep(0:2, 3)
out$region <- factor(out$region, levels = 0:2, labels = c("All respondents", "West Germany", "East Germany"))
out$treatment <- rep(1:3, each=3)
out$treatment <- factor(out$treatment, levels=1:3, labels = c("Control", "Economic insecurity", "Political underrepresentation"))

### define relevant statistics
# get means
out[1,3] <- mean(d$demsatis[d$treat == "control"], na.rm=T)
out[2,3] <- mean(d$demsatis[d$east == 0 & d$treat == "control"], na.rm=T)
out[3,3] <- mean(d$demsatis[d$east == 1 & d$treat == "control"], na.rm=T)
out[4,3] <- mean(d$demsatis[d$treat == "treat econ"], na.rm=T)
out[5,3] <- mean(d$demsatis[d$east == 0 & d$treat == "treat econ"], na.rm=T)
out[6,3] <- mean(d$demsatis[d$east == 1 & d$treat == "treat econ"], na.rm=T)
out[7,3] <- mean(d$demsatis[d$treat == "treat pol"], na.rm=T)
out[8,3] <- mean(d$demsatis[d$east == 0 & d$treat == "treat pol"], na.rm=T)
out[9,3] <- mean(d$demsatis[d$east == 1 & d$treat == "treat pol"], na.rm=T)

# estimate standard errors
out[1,4] <- sd(d$demsatis[d$treat == "control"], na.rm=T) / 
  sqrt(length(d$demsatis[d$treat == "control" & !is.na(d$demsatis) & !is.na(d$east)]))
out[2,4] <- sd(d$demsatis[d$east==0 & d$treat == "control"], na.rm=T) / 
  sqrt(length(d$demsatis[d$east==0 & d$treat == "control" & !is.na(d$demsatis) & !is.na(d$east)]))
out[3,4] <- sd(d$demsatis[d$east==1 & d$treat == "control"], na.rm=T) / 
  sqrt(length(d$demsatis[d$east==1 & d$treat == "control" & !is.na(d$demsatis) & !is.na(d$east)]))
out[4,4] <- sd(d$demsatis[d$treat == "treat econ"], na.rm=T) / 
  sqrt(length(d$demsatis[d$treat == "treat econ" & !is.na(d$demsatis) & !is.na(d$east)]))
out[5,4] <- sd(d$demsatis[d$east==0 & d$treat == "treat econ"], na.rm=T) / 
  sqrt(length(d$demsatis[d$east==0 & d$treat == "treat econ" & !is.na(d$demsatis) & !is.na(d$east)]))
out[6,4] <- sd(d$demsatis[d$east==1 & d$treat == "treat econ"], na.rm=T) / 
  sqrt(length(d$demsatis[d$east==1 & d$treat == "treat econ" & !is.na(d$demsatis) & !is.na(d$east)]))
out[7,4] <- sd(d$demsatis[d$treat == "treat pol"], na.rm=T) / 
  sqrt(length(d$demsatis[d$treat == "treat pol" & !is.na(d$demsatis) & !is.na(d$east)]))
out[8,4] <- sd(d$demsatis[d$east==0 & d$treat == "treat pol"], na.rm=T) / 
  sqrt(length(d$demsatis[d$east==0 & d$treat == "treat pol" & !is.na(d$demsatis) & !is.na(d$east)]))
out[9,4] <- sd(d$demsatis[d$east==1 & d$treat == "treat pol"], na.rm=T) / 
  sqrt(length(d$demsatis[d$east==1 & d$treat == "treat pol" & !is.na(d$demsatis) & !is.na(d$east)]))


### Plot
ggplot() + 
  geom_bar(data=out, aes(x=region, y=estimate - 1, fill=treatment), stat="identity", position=position_dodge(), width=.5) + 
  geom_errorbar(data=out, aes(x=region, ymin=estimate - 1 - 1.96 * se, ymax = estimate - 1 + 1.96*se, col=treatment),
                position=position_dodge(.5), width=.1, size=.75) + 
  scale_fill_manual(values = c("gray90", "gray70", "gray50"), name="") +
  scale_color_manual(values = c("gray60", "gray40", "gray20"), name="") + 
  scale_y_continuous(limits=c(0,3.1), breaks=seq(2,3.1,0.25),
                     labels = seq(2,3.1,0.25)+1) +
  coord_cartesian(ylim=c(2,3.1)) +
  theme_classic() + 
  xlab("") + 
  ylab("Democracy satisfaction\n(1=very dissatisfied; 6=very satisfied)") + 
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=12, face="bold"),
        axis.text.y = element_text(size=10),
        legend.text = element_text(size=12))





# Table 1 -----------------------------------------------------------------

### Regressions
m1 <- felm(demsatis ~ east + treat | 0 | 0 | ID, data=d)
m2 <- felm(demsatis ~ east*treat | 0 | 0 | ID, data=d)

### Stargazer output 
stargazer(m1,m2,
          omit.stat = c("rsq", "f", "ser"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(4),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),          
          covariate.labels=c("East Germany",
                             "Economic insecurity",
                             "Political underrepresentation",
                             "East Germany * Economic insecurity",
                             "East Germany * Political underrepresentation"), 
          digits = 3)





# Figure 3 ----------------------------------------------------------------

### code most important characteristic of democracy
d$dem_mostimportant <- 0
d$dem_mostimportant[d$dem_majoritarian_1 > d$dem_liberal_1 & 
                            d$dem_majoritarian_1 > d$dem_performance_1] <- 1
d$dem_mostimportant[d$dem_liberal_1 > d$dem_majoritarian_1 & 
                            d$dem_liberal_1 > d$dem_performance_1] <- 2
d$dem_mostimportant[d$dem_performance_1 > d$dem_majoritarian_1 & 
                            d$dem_performance_1 > d$dem_liberal_1] <- 3
d$dem_mostimportant <- factor(d$dem_mostimportant,
                                    levels = 0:3,
                                    labels = c("no clear preference", "majoritarian", "liberal", "performance"))

### create output container
out2 <- data.frame(matrix(nrow=26, ncol=5))
colnames(out2) <- c("class", "variable", "treat", "b", "se")
out2$class <- c(rep(1,10), rep(2,16))
out2$class <- factor(out2$class, levels=1:2, labels = c("Demographics", "Attitudes"))
out2$variable <- rep(1:13, each=2)
out2$variable <- factor(out2$variable, levels=1:13,
                        labels = c("Income quartiles", "Employment status", "Age quartiles", 
                                   "Education", "All demographics",
                                   "Economic satisfaction", "Political grievances", "Socialist attitudes",
                                   "Democracy: majoritarian", "Democracy: liberal", "Democracy: performance",
                                   "Democracy: most important", "All attitudes"))
out2$treat <- rep(1:2, 13)
out2$treat <- factor(out2$treat, levels=1:2, labels = c("Economic insecurity",
                                                        "Political underrepresentation"))

### Regressions
# adding baseline demographics
m1 <- felm(demsatis ~ east*treat + income_quartiles*treat | 0 | 0 | ID, data=d)
m2 <- felm(demsatis ~ east*treat + employment_status*treat | 0 | 0 | ID, data=d)
m3 <- felm(demsatis ~ east*treat + age_quartiles*treat | 0 | 0 | ID, data=d)
m4 <- felm(demsatis ~ east*treat + educ*treat | 0 | 0 | ID, data=d)
m5 <- felm(demsatis ~ east*treat + income_quartiles*treat + employment_status*treat + 
             age_quartiles*treat + educ*treat | 0 | 0 | ID, data=d)

# adding baseline attitudes
m6 <- felm(demsatis ~ east*treat + econsatis.today*treat | 0 | 0 | ID, data=d)
m7 <- felm(demsatis ~ east*treat + dontcare*treat | 0 | 0 | ID, data=d)
m8 <- felm(demsatis ~ east*treat + socialist.ideas*treat | 0 | 0 | ID, data=d)
m9 <- felm(demsatis ~ east*treat + dem_majoritarian_1*treat | 0 | 0 | ID, data=d)
m10 <- felm(demsatis ~ east*treat + dem_liberal_1*treat | 0 | 0 | ID, data=d)
m11 <- felm(demsatis ~ east*treat + dem_performance_1*treat | 0 | 0 | ID, data=d)
m12 <- felm(demsatis ~ east*treat + dem_mostimportant*treat | 0 | 0 | ID, data=d)
m13 <- felm(demsatis ~ east*treat + econsatis.today*treat + dontcare*treat + socialist.ideas*treat + 
              dem_majoritarian_1*treat + dem_liberal_1*treat + dem_performance_1*treat + dem_mostimportant*east | 0 | 0 | ID, data=d)


### extract coefficients
out2[1,4:5] <- summary(m1)$coef["east:treattreat econ",1:2]
out2[2,4:5] <- summary(m1)$coef["east:treattreat pol",1:2]
out2[3,4:5] <- summary(m2)$coef["east:treattreat econ",1:2]
out2[4,4:5] <- summary(m2)$coef["east:treattreat pol",1:2]
out2[5,4:5] <- summary(m3)$coef["east:treattreat econ",1:2]
out2[6,4:5] <- summary(m3)$coef["east:treattreat pol",1:2]
out2[7,4:5] <- summary(m4)$coef["east:treattreat econ",1:2]
out2[8,4:5] <- summary(m4)$coef["east:treattreat pol",1:2]
out2[9,4:5] <- summary(m5)$coef["east:treattreat econ",1:2]
out2[10,4:5] <- summary(m5)$coef["east:treattreat pol",1:2]
out2[11,4:5] <- summary(m6)$coef["east:treattreat econ",1:2]
out2[12,4:5] <- summary(m6)$coef["east:treattreat pol",1:2]
out2[13,4:5] <- summary(m7)$coef["east:treattreat econ",1:2]
out2[14,4:5] <- summary(m7)$coef["east:treattreat pol",1:2]
out2[15,4:5] <- summary(m8)$coef["east:treattreat econ",1:2]
out2[16,4:5] <- summary(m8)$coef["east:treattreat pol",1:2]
out2[17,4:5] <- summary(m9)$coef["east:treattreat econ",1:2]
out2[18,4:5] <- summary(m9)$coef["east:treattreat pol",1:2]
out2[19,4:5] <- summary(m10)$coef["east:treattreat econ",1:2]
out2[20,4:5] <- summary(m10)$coef["east:treattreat pol",1:2]
out2[21,4:5] <- summary(m11)$coef["east:treattreat econ",1:2]
out2[22,4:5] <- summary(m11)$coef["east:treattreat pol",1:2]
out2[23,4:5] <- summary(m12)$coef["east:treattreat econ",1:2]
out2[24,4:5] <- summary(m12)$coef["east:treattreat pol",1:2]
out2[25,4:5] <- summary(m13)$coef["east:treattreat econ",1:2]
out2[26,4:5] <- summary(m13)$coef["east:treattreat pol",1:2]

### Plot
ggplot() +
  geom_hline(yintercept = 0, lty = 2) +
  geom_point(data = out2, aes(x=variable, y=b, shape=treat, col=treat),
             position = position_dodge(width=0.8), size=2.2) + 
  geom_errorbar(data = out2, aes(x=variable, ymin = b - 1.96*se, ymax = b+1.96*se, col=treat), width=0,
                position = position_dodge(width=0.8), size=0.8) +
  scale_color_manual(values = c("gray10", "gray70")) +
  scale_x_discrete(limits=rev(levels(out2$variable))) +
  coord_flip() + 
  #facet_grid(class~treat, scales="free_y") + 
  xlab("Included: Interaction between\n treatment groups and...") +
  ylab("Estimated difference in the treatment effect\nbetween West and East Germans") +
  theme_classic() +
  theme(axis.title = element_text(size=12),
        axis.text.y = element_text(size=11),
        axis.text.x = element_text(size=9),
        strip.text = element_text(size=13, face="bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size=12))




# Table 2 -----------------------------------------------------------------

### manually obtain values for each cell
out <- matrix(nrow=3, ncol=8)
out[1,1] <- length(d$treat_applicable[d$east==0 & d$treat_applicable=="Control"])
out[1,4] <- length(d$treat_applicable[d$east==1 & d$treat_applicable=="Control"])
out[2,1] <- length(d$treat_applicable[d$east==0 & d$treat=="treat econ"])
out[2,2] <- length(d$treat_applicable[d$east==0 & d$treat_applicable=="Econ yes"])
out[2,3] <- round((out[2,2] / out[2,1])*100, digits=2)
out[2,4] <- table(d$east==1 & d$treat=="treat econ")[2]
out[2,5] <- length(d$treat_applicable[d$east==1 & d$treat_applicable=="Econ yes"])
out[2,6] <- round((out[2,5] / out[2,4])*100, digits=2)
out[2,7] <- round(out[2,3] - out[2,6], digits=2)
out[2,8] <- round((t.test(as.numeric(d$treat_applicable[d$treat == "treat econ" & d$east == 0]),
                          as.numeric(d$treat_applicable[d$treat == "treat econ" & d$east == 1]))$p.value), digits=3)

out[3,1] <- length(d$treat_applicable[d$east==0 & d$treat=="treat pol"])
out[3,2] <- length(d$treat_applicable[d$east==0 & d$treat_applicable=="Pol yes"])
out[3,3] <- round((out[3,2] / out[3,1])*100, digits=2)
out[3,4] <- table(d$east==1 & d$treat=="treat pol")[2]
out[3,5] <- length(d$treat_applicable[d$east==1 & d$treat_applicable=="Pol yes"])
out[3,6] <- round((out[3,5] / out[3,4])*100, digits=2)
out[3,7] <- round(out[3,3] - out[3,6], digits=2)
out[3,8] <- round((t.test(as.numeric(d$treat_applicable[d$treat == "treat pol" & d$east == 0]),
                          as.numeric(d$treat_applicable[d$treat == "treat pol" & d$east == 1]))$p.value), digits=3)

### define column and row names
rownames(out) <- c("Control", "Political", "Economic")
colnames(out) <- c("East Total", "East applicable N", "East applicable share", 
                   "West Total", "West applicable N"," West applicable share",
                   "difference ppt", "difference p-value")


### Table output
print(xtable(out, caption = "East and West Germans are equally likely to have experienced insecurity or under-representation"),
      caption.placement = "top")




# Table 3 -----------------------------------------------------------------

## Preparing the word list
# load word list
wordlist <- read_xlsx("Data/SentiWS_v2/SentiWS_v2.0_All.xlsx")

# all words lower case
wordlist$term <- tolower(wordlist$term)

# stemming the words
wordlist.stemmed <- tm_map(Corpus(VectorSource(wordlist$term)), stemDocument, language = "german")

# putting the words back into the dataframes
wordlist$stemmed <- data.frame(text = get("content", wordlist.stemmed))
wordlist$stemmed <- as.matrix(wordlist$stemmed)


## Preparing response text for analysis
# create the responses
d$response <- NA
d$response[d$treat_applicable == "Econ yes"] <- 
  as.character(d$treat_econ[d$treat_applicable == "Econ yes"])
d$response[d$treat_applicable == "Pol yes"] <- 
  as.character(d$treat_pol[d$treat_applicable == "Pol yes"])

# Stemming the responses
aux <- iconv(d$response, from="UTF-8", to="ASCII//TRANSLIT", sub="byte")

# change Umlaute
aux <- gsub('\\"u', 'ü', aux)
aux <- gsub('\\"a', 'ä', aux)
aux <- gsub('\\"o', 'ö', aux)

# convert tweets into corpus
docs.germany <- Corpus(VectorSource(aux))

# move all words to lower case
docs.germany <- tm_map(docs.germany,content_transformer(tolower))

# remove potentially problematic symbols
toSpace <- content_transformer(function(x, pattern) { return (gsub(pattern, " ", x))})
docs.germany <- tm_map(docs.germany, toSpace, "-")
docs.germany <- tm_map(docs.germany, toSpace, "'")
docs.germany <- tm_map(docs.germany, toSpace, "`")
docs.germany <- tm_map(docs.germany, toSpace, "•")
docs.germany <- tm_map(docs.germany, toSpace, "``")
docs.germany <- tm_map(docs.germany, toSpace, "''")
docs.germany <- tm_map(docs.germany, toSpace, "\n")
docs.germany <- tm_map(docs.germany, toSpace, ">")
docs.germany <- tm_map(docs.germany, toSpace, "\\.")
docs.germany <- tm_map(docs.germany, toSpace, "\\,")
docs.germany <- tm_map(docs.germany, toSpace, "amp")

# remove punctuation and stop words
docs.germany <- tm_map(docs.germany, removePunctuation)
docs.germany <- tm_map(docs.germany, removeNumbers)
docs.germany <- tm_map(docs.germany, removeWords, stopwords("german"))
myStopwords <- c("dass", "d", "g", "u", "s")
docs.germany <- tm_map(docs.germany, removeWords, myStopwords)
docs.germany <- tm_map(docs.germany, stripWhitespace)

# word stemming
d$response.stemmed <- stemDocument(docs.germany$content, language="german")

# code length of response
d$response.length <- str_count(d$response, "\\S+")


### compute the sentiment score for each response
d$sentiment1 <- NA
d$sentiment2 <- NA
d[,c("sentiment1", "sentiment2")] <- t(apply(as.matrix(d$response.stemmed), MARGIN=1, FUN = function(x){
  word.list <- str_split(x, ' ')
  words <- unlist(word.list)
  matches <- match(words, wordlist$stemmed)
  score <- wordlist$score[matches]
  mean1 <- mean(score, na.rm=T)
  mean1[is.nan(mean1)] <- 0
  score[is.na(score)] <- 0
  mean2 <- mean(score)
  return(c(mean1, mean2))
}))



### create dichotomous variables measuring the content of responses
# Econ: income
d$treat_econ_type_income <- NA
d$treat_econ_type_income[d$treat_econ_type == "income"] <- 1
d$treat_econ_type_income[d$treat_econ_type != "income" & d$treat == "treat econ" & d$treat_applicable == "Econ yes"] <- 0

# Econ: unemployed
d$treat_econ_type_unemp <- NA
d$treat_econ_type_unemp[d$treat_econ_type == "unemployed"] <- 1
d$treat_econ_type_unemp[d$treat_econ_type != "unemployed" & d$treat == "treat econ" & d$treat_applicable == "Econ yes"] <- 0

# Econ: feeling
d$treat_econ_type_feeling <- NA
d$treat_econ_type_feeling[d$treat_econ_type == "feeling"] <- 1
d$treat_econ_type_feeling[d$treat_econ_type != "feeling" & d$treat == "treat econ" & d$treat_applicable == "Econ yes"] <- 0


# Pol: economy
d$treat_pol_type_econ <- NA
d$treat_pol_type_econ[d$treat_pol_type == "economy"] <- 1
d$treat_pol_type_econ[d$treat_pol_type != "economy" & d$treat == "treat pol" & d$treat_applicable == "Pol yes"] <- 0

# Pol: policy differences
d$treat_pol_type_policy <- NA
d$treat_pol_type_policy[d$treat_pol_type == "policy differences"] <- 1
d$treat_pol_type_policy[d$treat_pol_type != "policy differences" & d$treat == "treat pol" & d$treat_applicable == "Pol yes"] <- 0

# Pol: representation
d$treat_pol_type_rep <- NA
d$treat_pol_type_rep[d$treat_pol_type == "representation"] <- 1
d$treat_pol_type_rep[d$treat_pol_type != "representation" & d$treat == "treat pol" & d$treat_applicable == "Pol yes"] <- 0

# Pol: feeling
d$treat_pol_type_feeling <- NA
d$treat_pol_type_feeling[d$treat_pol_type == "feeling"] <- 1
d$treat_pol_type_feeling[d$treat_pol_type != "feeling" & d$treat == "treat pol" & d$treat_applicable == "Pol yes"] <- 0


### create output container
# empty container
out <- matrix(nrow=11, ncol=4)

# Get distribution of topics
out[1,1] <- round(mean(d$treat_econ_type_income[d$east==0], na.rm=T)*100, digits=2)
out[1,2] <- round(mean(d$treat_econ_type_income[d$east==1], na.rm=T)*100, digits=2)
out[1,3] <- out[1,2] - out[1,1]
out[1,4] <- round(coef(summary(felm(treat_econ_type_income ~ east | 0 | 0 | ID, data=d)))[2,4], digits=2)
out[2,1] <- round(mean(d$treat_econ_type_unemp[d$east==0], na.rm=T)*100, digits=2)
out[2,2] <- round(mean(d$treat_econ_type_unemp[d$east==1], na.rm=T)*100, digits=2)
out[2,3] <- out[2,2] - out[2,1]
out[2,4] <- round(coef(summary(felm(treat_econ_type_unemp ~ east | 0 | 0 | ID, data=d)))[2,4], digits=2)
out[3,1] <- round(mean(d$treat_econ_type_feeling[d$east==0], na.rm=T)*100, digits=2)
out[3,2] <- round(mean(d$treat_econ_type_feeling[d$east==1], na.rm=T)*100, digits=2)
out[3,3] <- out[3,2] - out[3,1]
out[3,4] <- round(coef(summary(felm(treat_econ_type_feeling ~ east | 0 | 0 | ID, data=d)))[2,4], digits=2)
out[4,1] <- round(mean(d$treat_pol_type_econ[d$east==0], na.rm=T)*100, digits=2)
out[4,2] <- round(mean(d$treat_pol_type_econ[d$east==1], na.rm=T)*100, digits=2)
out[4,3] <- out[4,2] - out[4,1]
out[4,4] <- round(coef(summary(felm(treat_pol_type_econ ~ east | 0 | 0 | ID, data=d)))[2,4], 2)
out[5,1] <- round(mean(d$treat_pol_type_policy[d$east==0], na.rm=T)*100, digits=2)
out[5,2] <- round(mean(d$treat_pol_type_policy[d$east==1], na.rm=T)*100, digits=2)
out[5,3] <- out[5,2] - out[5,1]
out[5,4] <- round(coef(summary(felm(treat_pol_type_policy ~ east | 0 | 0 | ID, data=d)))[2,4],digits=2)
out[6,1] <- round(mean(d$treat_pol_type_rep[d$east==0], na.rm=T)*100, digits=2)
out[6,2] <- round(mean(d$treat_pol_type_rep[d$east==1], na.rm=T)*100, digits=2)
out[6,3] <- out[6,2] - out[6,1]
out[6,4] <- round(coef(summary(felm(treat_pol_type_rep ~ east | 0 | 0 | ID, data=d)))[2,4],digits=2)
out[7,1] <- round(mean(d$treat_pol_type_feeling[d$east==0], na.rm=T)*100, digits=2)
out[7,2] <- round(mean(d$treat_pol_type_feeling[d$east==1], na.rm=T)*100, digits=2)
out[7,3] <- out[7,2] - out[7,1]
out[7,4] <- round(coef(summary(felm(treat_pol_type_feeling ~ east | 0 | 0 | ID, data=d)))[2,4], digits=2)


# Add number of characters
out[8,1] <- round(mean(nchar(as.character(d$treat_econ[d$east==0 & d$treat_applicable == "Econ yes"])), na.rm=T), digits=2)
out[8,2] <- round(mean(nchar(as.character(d$treat_econ[d$east==1 & d$treat_applicable == "Econ yes"])), na.rm=T), digits=2)
out[8,3] <- out[8,2] - out[8,1]
out[8,4] <- round(coef(summary(felm(nchar(as.character(treat_econ)) ~ east | 0 | 0 | ID, data=d, subset=treat_applicable == "Econ yes")))[2,4], digits=2)

out[9,1] <- round(mean(nchar(as.character(d$treat_pol[d$east==0 & d$treat_applicable == "Pol yes"])), na.rm=T), digits=2)
out[9,2] <- round(mean(nchar(as.character(d$treat_pol[d$east==1 & d$treat_applicable == "Pol yes"])), na.rm=T), digits=2)
out[9,3] <- out[9,2] - out[9,1]
out[9,4] <- round(coef(summary(felm(nchar(as.character(treat_pol)) ~ east | 0 | 0 | ID, data=d, subset=treat_applicable == "Pol yes")))[2,4], digits=2)


# Sentiment analysis
out[10,1] <- round(mean(d$sentiment1[d$treat_applicable == "Econ yes" & d$east==0], na.rm=T), digits=2)
out[10,2] <- round(mean(d$sentiment1[d$treat_applicable == "Econ yes" & d$east==1], na.rm=T), digits=2)
out[10,3] <- out[10,2] - out[10,1]
out[10,4] <- round(coef(summary(felm(sentiment1 ~ east | 0 | 0 | ID, data=d, subset=treat_applicable == "Econ yes")))[2,4], digits=2)
out[11,1] <- round(mean(d$sentiment1[d$treat_applicable == "Pol yes" & d$east==0], na.rm=T), digits=2)
out[11,2] <- round(mean(d$sentiment1[d$treat_applicable == "Pol yes" & d$east==1], na.rm=T), digits=2)
out[11,3] <- out[11,2] - out[11,1]
out[11,4] <- round(coef(summary(felm(sentiment1 ~ east | 0 | 0 | ID, data=d, subset=treat_applicable == "Pol yes")))[2,4], digits=2)



# Turn into data frame
out <- data.frame(out)
colnames(out) <- c("West", "East", "Diff", "pvalue")
out$variable <- c("Economic: Income", "Economic: Unemployment", "Economic: Feeling",
                  "Political: Economic", "Political: Policy differences", "Political: Representation", "Political: Feeling",
                  "Economic insecurity", "Political under-representation",
                  "Economic insecurity", "Political under-representation")
out <- out[,c("variable", "West", "East", "Diff","pvalue")]


### output table
print(xtable(out, caption = "East and West Germans report qualitatively similar experiences"),
      caption.placement = "top", include.rownames=FALSE)





# Table 4 -----------------------------------------------------------------

### Regresssions
m1 <- felm(demsatis ~ treat_applicable | 0 | 0 | ID, data=d)
m2 <- felm(demsatis ~ treat_applicable | 0 | 0 | ID, data=d, subset=east==0)
m3 <- felm(demsatis ~ treat_applicable | 0 | 0 | ID, data=d, subset=east==1)

### Stargazer output 
stargazer(m1, m2, m3,
          omit.stat = c("rsq", "f", "ser"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(3),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),          
          covariate.labels=c("Economic insecurity -- no relevant experience",
                             "Economic insecurity -- relevant experience",
                             "Political underrepresentation -- no relevant experience",
                             "Political underrepresentation -- relevant experience"), 
          digits = 3)





# Figure A1 ---------------------------------------------------------------

### load and prepare shapefile for Kreise
# load
shp1 <- readOGR("Data/Shapefile Kreis 2018/VG1000_F.shp")

# select the correct units
shp1 <- subset(shp1, shp1@data$GF==4)

# change CRS
shp1 <- spTransform(shp1, CRS("+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"))

# fortify
shp.kreise <- fortify(shp1, region="RS")

### load and prepare shapefile for Land
# load
shp2 <- readOGR("Data/Shapefile Land 2018/VG1000_LAN.shp")

# select the correct units
shp2 <- subset(shp2, shp2@data$GF==4)

# fortify
shp.land <- fortify(shp2, region="AGS")

### adjust lon/lat to map CRS
# extract lon/lat
sp <- SpatialPoints(d[!is.na(d$lonNEW),c("lonNEW", "latNEW")], CRS("+proj=longlat +datum=WGS84"))

# change CRS
sp2 <- spTransform(sp, CRS("+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"))

# turn into data frame
sp2 <- data.frame(sp2)

### map
ggplot() + 
  geom_polygon(data=shp.kreise, aes(x=long, y=lat, group=group), col="gray80", fill="white", size=.15) + 
  geom_path(data=shp.land, aes(x=long, y=lat, group=group), col="gray65", size=.5) + 
  geom_jitter(data=sp2, aes(x=lonNEW, y=latNEW), size=.8, alpha=.7, width=10000, height=10000, col="salmon4") +
  coord_cartesian(xlim=c(250000,950000), ylim=c(5200000, 6100000)) +
  theme_bw() + 
  theme(axis.title=element_blank(), 
        axis.text=element_blank(), 
        axis.ticks=element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.title = element_text(size=20, face="bold", hjust=.5),
        panel.background = element_rect(fill = 'gray97', colour = NA),
        legend.position="bottom",
        legend.text = element_text(size=16, face = "bold"))  




# Figure A2 ---------------------------------------------------------------


##### Definition of democracy
### create averages
# output container
out <- matrix(nrow = 6, ncol=4)
out[,4] <- c("West", "East", "West","East","West","East")
out[,3] <- c(1,1,2,2,3,3)

# Majoritarian democracy
out[1,1] <- mean(d$dem_majoritarian_1[d$east==0], na.rm=T)
out[2,1] <- mean(d$dem_majoritarian_1[d$east==1], na.rm=T)
out[1,2] <- sd(d$dem_majoritarian_1[d$east==0], na.rm=T) / 
  sqrt(length(d$dem_majoritarian_1[d$east==0 & !is.na(d$dem_majoritarian_1 & !is.na(d$east))]))
out[2,2] <- sd(d$dem_majoritarian_1[d$east==1], na.rm=T) / 
  sqrt(length(d$dem_majoritarian_1[d$east==1 & !is.na(d$dem_majoritarian_1 & !is.na(d$east))]))

# Liberal democracy
out[3,1] <- mean(d$dem_liberal_1[d$east==0], na.rm=T)
out[4,1] <- mean(d$dem_liberal_1[d$east==1], na.rm=T)
out[3,2] <- sd(d$dem_liberal_1[d$east==0], na.rm=T) / 
  sqrt(length(d$dem_liberal_1[d$east==0 & !is.na(d$dem_liberal_1 & !is.na(d$east))]))
out[4,2] <- sd(d$dem_liberal_1[d$east==1], na.rm=T) / 
  sqrt(length(d$dem_liberal_1[d$east==1 & !is.na(d$dem_liberal_1 & !is.na(d$east))]))

# Performance democracy
out[5,1] <- mean(d$dem_performance_1[d$east==0], na.rm=T)
out[6,1] <- mean(d$dem_performance_1[d$east==1], na.rm=T)
out[5,2] <- sd(d$dem_performance_1[d$east==0], na.rm=T) / 
  sqrt(length(d$dem_performance_1[d$east==0 & !is.na(d$dem_performance_1 & !is.na(d$east))]))
out[6,2] <- sd(d$dem_performance_1[d$east==1], na.rm=T) / 
  sqrt(length(d$dem_performance_1[d$east==1 & !is.na(d$dem_performance_1 & !is.na(d$east))]))

# convert to data frame
out <- data.frame(out)
colnames(out) <- c("mean", "se", "type", "region")
out$mean <- as.numeric(as.character(out$mean))
out$se <- as.numeric(as.character(out$se))
out$type <- factor(out$type, levels=c(1,2,3),
                   labels = c("Majoritarian", "Liberal", "Performance"))


### Plot
ggplot() + 
  geom_bar(data=out, aes(x=type, y=mean, fill=region), stat="identity", position=position_dodge(), width=.5) + 
  geom_errorbar(data=out, aes(x=type, ymin=mean - 1.96 * se, ymax = mean + 1.96*se, col=region),
                position=position_dodge(.5), width=.1, size=.75) + 
  scale_fill_manual(values = c("gray90", "gray70"), name="") +
  scale_color_manual(values = c("gray60", "gray40"), name="") + 
  scale_y_continuous(limits=c(0,10), breaks=seq(0,10,2)) +
  theme_classic() + 
  xlab("Definition of Democracy") + 
  ylab("Importance\n(0=not at all important; 10=essential)") + 
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=12, face="bold"),
        axis.text.y = element_text(size=10),
        legend.text = element_text(size=12))




##### Attitudes
### create averages
# output container
out <- matrix(nrow = 6, ncol=4)
out[,4] <- c("West", "East", "West","East","West","East")
out[,3] <- c(1,1,2,2,3,3)

# Political representation
out[1,1] <- mean(d$dontcare[d$east==0], na.rm=T)
out[2,1] <- mean(d$dontcare[d$east==1], na.rm=T)
out[1,2] <- sd(d$dontcare[d$east==0], na.rm=T) / 
  sqrt(length(d$dontcare[d$east==0 & !is.na(d$dontcare & !is.na(d$east))]))
out[2,2] <- sd(d$dontcare[d$east==1], na.rm=T) / 
  sqrt(length(d$dontcare[d$east==1 & !is.na(d$dontcare & !is.na(d$east))]))

# Economic satisfaction
out[3,1] <- mean(d$econsatis.today[d$east==0], na.rm=T)
out[4,1] <- mean(d$econsatis.today[d$east==1], na.rm=T)
out[3,2] <- sd(d$econsatis.today[d$east==0], na.rm=T) / 
  sqrt(length(d$econsatis.today[d$east==0 & !is.na(d$econsatis.today & !is.na(d$east))]))
out[4,2] <- sd(d$econsatis.today[d$east==1], na.rm=T) / 
  sqrt(length(d$econsatis.today[d$east==1 & !is.na(d$econsatis.today & !is.na(d$east))]))

# Socialist legacies
out[5,1] <- mean(d$socialist.ideas[d$east==0], na.rm=T)
out[6,1] <- mean(d$socialist.ideas[d$east==1], na.rm=T)
out[5,2] <- sd(d$socialist.ideas[d$east==0], na.rm=T) / 
  sqrt(length(d$socialist.ideas[d$east==0 & !is.na(d$socialist.ideas & !is.na(d$east))]))
out[6,2] <- sd(d$socialist.ideas[d$east==1], na.rm=T) / 
  sqrt(length(d$socialist.ideas[d$east==1 & !is.na(d$socialist.ideas & !is.na(d$east))]))

# convert to data frame
out <- data.frame(out)
colnames(out) <- c("mean", "se", "type", "region")
out$type <- factor(out$type, levels=c(1,2,3),
                   labels = c("Politicians\ndon't care", "Economic\nsatisfaction", "Socialist\nideas"))
out$mean <- as.numeric(as.character(out$mean))
out$se <- as.numeric(as.character(out$se))


### Plot
ggplot() + 
  geom_bar(data=out, aes(x=type, y=mean, fill=region), stat="identity", position=position_dodge(), width=.5) + 
  geom_errorbar(data=out, aes(x=type, ymin=mean - 1.96 * se, ymax = mean + 1.96*se, col=region),
                position=position_dodge(.5), width=.1, size=.75) + 
  scale_fill_manual(values = c("gray90", "gray70"), name="") +
  scale_color_manual(values = c("gray60", "gray40"), name="") + 
  scale_y_continuous(limits=c(0,1.2), breaks=seq(0,1.2,0.2)) +
  theme_classic() + 
  xlab("Definition of Democracy") + 
  ylab("Agreement\n(-2=strong disagreement; 2=strong agreement)") + 
  theme(legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size=14, face="bold"),
        axis.text.x = element_text(size=12, face="bold"),
        axis.text.y = element_text(size=10),
        legend.text = element_text(size=12))





# Table A1 ----------------------------------------------------------------

### Write a function that returns summary statistics
dFun <- function(dv){
  # get values
  mean <- mean(d[,dv], na.rm=T)
  sd <- sd(d[,dv], na.rm=T)
  n <- length(d[!is.na(d[,dv]),dv])
  # put together
  out <- c(mean, sd, n)
  # output
  return(out)
}

### create table
descriptives <- rbind(dFun(dv="demsatis"), 
                      dFun(dv="female"),
                      dFun(dv="age2"),
                      dFun(dv="east"),
                      dFun(dv="income.numeric"),
                      dFun(dv="educ_89"),
                      dFun(dv="educ_10"),
                      dFun(dv="educ_12"),
                      dFun(dv="employment_work"),
                      dFun(dv="employment_notwork"),
                      dFun(dv="employment_homemaker"),
                      dFun(dv="employment_inschool"),
                      dFun(dv="employment_retired"),
                      dFun(dv="econsatis.today"),
                      dFun(dv="dontcare"),
                      dFun(dv="socialist.ideas"),
                      dFun(dv="dem_majoritarian_1"),
                      dFun(dv="dem_liberal_1"),
                      dFun(dv="dem_performance_1"))
colnames(descriptives) <- c("Mean", "SD", "N")
rownames(descriptives) <- c("Democracy satisfaction", "Female","Age","East", "Income", 
                            "Education: 8-9 years", "Education: 10 years", "Education: 12/13 years", 
                            "Employment: employed", "Employment: not working", 
                            "Employment: homemaker", "Employment: in school", "Employment: retired",
                            "Economic satisfaction", "Political grievances", "Socialist attitudes",
                            "Democracy: majoritarian", "Democracy: liberal", "Democracy: performance")

### Table output
print(xtable(descriptives, digits=2, caption = "Descriptive statistics"), 
      caption.placement = "top")




# Table A2 ----------------------------------------------------------------

### get relevant statistics
xtable(table(d$treat[!is.na(d$demsatis)], d$east[!is.na(d$demsatis)]), digits=2)
xtable(prop.table(table(d$treat[!is.na(d$demsatis)], d$east[!is.na(d$demsatis)]))*100, digits=2)




# Table A3 ----------------------------------------------------------------

### Write a function that runs a balance test
balance <- function(dv){
  # run regression
  reg <- felm(d[,dv] ~ as.factor(treat) | 0 | 0 | ID, data=d)
  # extract values
  b1 <- summary(reg)$coef[2,1]
  p1 <- summary(reg)$coef[2,4]
  b2 <- summary(reg)$coef[3,1]
  p2 <- summary(reg)$coef[3,4]
  
  # get means
  mean1 <- mean(d[d$treat == "control",dv], na.rm=T)
  mean2 <- mean(d[d$treat == "treat econ",dv], na.rm=T)
  mean3 <- mean(d[d$treat == "treat pol",dv], na.rm=T)
  # put together
  out <- c(mean1, mean2, mean3, b1,p1,b2,p2)
  # output
  return(out)
}

balancetable <- rbind(balance(dv="female"),
                      balance(dv="age2"),
                      balance(dv="east"),
                      balance(dv="income.numeric"),
                      balance(dv="educ_89"),
                      balance(dv="educ_10"),
                      balance(dv="educ_12"),
                      balance(dv="employment_work"),
                      balance(dv="employment_notwork"),
                      balance(dv="employment_homemaker"),
                      balance(dv="employment_inschool"),
                      balance(dv="employment_retired"),
                      balance(dv="econsatis.today"),
                      balance(dv="dontcare"),
                      balance(dv="socialist.ideas"),
                      balance(dv="dem_majoritarian_1"),
                      balance(dv="dem_liberal_1"),
                      balance(dv="dem_performance_1"))
colnames(balancetable) <- c("Mean Control", "Mean Econ", "Mean Treat Pol", "Econ Diff", "Econ p", "Pol Diff", "Pol p")
rownames(balancetable) <- c("Female","Age","East", "Income", "Education: 8-9 years", "Education: 10 years",
                            "Education: 12/13 years", "Employment: employed", "Employment: not working", 
                            "Employment: homemaker", "Employment: in school", "Employment: retired",
                            "Economic satisfaction", "Political grievances", "Socialist attitudes",
                            "Democracy: majoritarian", "Democracy: liberal", "Democracy: performance")


### Table output
print(xtable(balancetable, digits=2, caption = "Balance statistics"),
      caption.placement = "top")





# Table A4 ----------------------------------------------------------------

### Regressions
m1 <- felm(demsatis ~ east*treat + income_quartiles*treat | 0 | 0 | ID, data=d)
m2 <- felm(demsatis ~ east*treat + employment_status*treat | 0 | 0 | ID, data=d)
m3 <- felm(demsatis ~ east*treat + age_quartiles*treat | 0 | 0 | ID, data=d)
m4 <- felm(demsatis ~ east*treat + educ*treat | 0 | 0 | ID, data=d)
m5 <- felm(demsatis ~ east*treat + income_quartiles*treat + employment_status*treat + 
             age_quartiles*treat + educ*treat | 0 | 0 | ID, data=d)

### Stargazer output
stargazer(m1, m2, m3, m4, m5,
          omit.stat = c("rsq", "f", "ser"),
          omit = c("income_quartiles", "employment_status", "age_quartiles", "educ"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(6),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),          
          covariate.labels=c("East Germany",
                             "Economic insecurity",
                             "Political underrepresentation",
                             "East Germany * Economic insecurity",
                             "East Germany * Political underrepresentation"), 
          digits = 3,
          add.lines = list(
            c("Income quartile * Treatments?", "YES", "NO", "NO", "NO", "YES"),
            c("Employment status * Treatments?", "NO", "YES", "NO", "NO", "YES"),
            c("Age quartile * Treatments?", "NO", "NO", "YES", "NO", "YES"),
            c("Education * Treatments?", "NO", "NO", "NO", "YES", "YES")))






# Table A5 ----------------------------------------------------------------

### Regressions
m6 <- felm(demsatis ~ east*treat + econsatis.today*treat | 0 | 0 | ID, data=d)
m7 <- felm(demsatis ~ east*treat + dontcare*treat | 0 | 0 | ID, data=d)
m8 <- felm(demsatis ~ east*treat + socialist.ideas*treat | 0 | 0 | ID, data=d)
m9 <- felm(demsatis ~ east*treat + dem_majoritarian_1*treat | 0 | 0 | ID, data=d)
m10 <- felm(demsatis ~ east*treat + dem_liberal_1*treat | 0 | 0 | ID, data=d)
m11 <- felm(demsatis ~ east*treat + dem_performance_1*treat | 0 | 0 | ID, data=d)
m12 <- felm(demsatis ~ east*treat + dem_mostimportant*treat | 0 | 0 | ID, data=d)
m13 <- felm(demsatis ~ east*treat + econsatis.today*treat + dontcare*treat + socialist.ideas*treat + 
              dem_majoritarian_1*treat + dem_liberal_1*treat + dem_performance_1*treat + dem_mostimportant*east | 0 | 0 | ID, data=d)

### Stargazer output
stargazer(m6,m7,m8,m9,m10,m11,m12,m13,
          omit.stat = c("rsq", "f", "ser"),
          omit = c("econsatis.today", "dontcare", "socialist.ideas", "dem_majoritarian_1", 
                   "dem_liberal_1", "dem_performance_1", "dem_mostimportant"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(8),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),          
          covariate.labels=c("East Germany",
                             "Economic insecurity",
                             "Political underrepresentation",
                             "East Germany * Economic insecurity",
                             "East Germany * Political underrepresentation"), 
          digits = 3,
          add.lines = list(
            c("Economic satisfaction * Treatments?", "YES", "NO", "NO", "NO", "NO", "NO", "NO", "YES"),
            c("Political grievances * Treatments?", "NO", "YES", "NO", "NO", "NO", "NO", "NO","YES"),
            c("Socialist attitudes * Treatments?", "NO", "NO", "YES", "NO", "NO", "NO", "NO", "YES"),
            c("Democracy: majoritarian * Treatments?", "NO", "NO", "NO", "YES", "NO", "NO", "NO", "YES"),
            c("Democracy: liberal * Treatments?", "NO", "NO", "NO", "NO", "YES", "NO", "NO", "YES"),
            c("Democracy: performance * Treatments?", "NO", "NO", "NO", "NO", "NO", "YES", "NO", "YES"),
            c("Democracy: most important * Treatments?", "NO", "NO", "NO", "NO", "NO", "NO", "YES", "YES")))





# Table A6 ----------------------------------------------------------------

### coding content
d$content <- NA
d$content[d$treat == "control"] <- 0
d$content[d$treat_econ_type_income == 1] <- 1
d$content[d$treat_econ_type_unemp == 1] <- 2
d$content[d$treat_econ_type_feeling == 1] <- 3
d$content[d$treat_pol_type_econ == 1] <- 4
d$content[d$treat_pol_type_policy == 1] <- 5
d$content[d$treat_pol_type_rep == 1] <- 6
d$content[d$treat_pol_type_feeling == 1] <- 7
d$content <- factor(d$content, levels=0:7,
                          labels = c("control", "econ: income", "econ: unemployment", "econ: feeling",
                                     "pol: economic", "pol: policy", "pol: representation", "pol: feeling"))

### coding length
d$length <- NA
d$length[d$treat == "control"] <- 0
d$length[d$treat != "control"] <- d$response.length[d$treat != "control"]

### coding sentiment
d$sentimentAUX <- NA
d$sentimentAUX[d$treat == "control"] <- 0
d$sentimentAUX[d$treat !="control"] <- d$sentiment1[d$treat !="control"]


### Regressions
m1 <- felm(demsatis ~ east * treat + content | 0 | 0 | ID, data=d, subset=treat_applicable %in% c("Control", "Econ yes", "Pol yes"))
m2 <- felm(demsatis ~ east * treat + length | 0 | 0 | ID, data=d, subset=treat_applicable %in% c("Control", "Econ yes", "Pol yes"))
m3 <- felm(demsatis ~ east * treat + sentimentAUX | 0 | 0 | ID, data=d, subset=treat_applicable %in% c("Control", "Econ yes", "Pol yes"))


### Stargazer output 
stargazer(m1, m2, m3,
          omit.stat = c("rsq", "f", "ser"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(4),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),       
          covariate.labels=c("East Germany",
                             "Economic insecurity",
                             "Political underrepresentation",
                             "Content: econ--income",
                             "Content: econ--unemployment",
                             "Content: econ--feeling",
                             "Content: pol--economic",
                             "Content: pol--policy",
                             "Content: pol--representation",
                             "Content: pol--feeling",
                             "Length",
                             "Sentiment score",
                             "East Germany * Economic insecurity",
                             "East Germany * Political underrepresentation"), 
          digits = 3)




# Table A7 ----------------------------------------------------------------

### Regressions
m1 <- felm(demsatis ~ east + treat | 0 | 0 | ID, data=d, subset=stateNEW != "Berlin")
m2 <- felm(demsatis ~ east*treat | 0 | 0 | ID, data=d, subset=stateNEW != "Berlin")


### Stargazer output 
stargazer(m1, m2,
          omit.stat = c("rsq", "f", "ser"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(2),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),          
          covariate.labels=c("East Germany",
                             "Economic insecurity",
                             "Political underrepresentation",
                             "East Germany * Economic insecurity",
                             "East Germany * Political underrepresentation"), 
          digits = 3)




# Table A8 ----------------------------------------------------------------

### Regressions
m1 <- polr(as.factor(demsatis) ~ east + treat, data = d, Hess=TRUE, method="logistic")
m2 <- polr(as.factor(demsatis) ~ east + treat, data = d, Hess=TRUE, method="probit")

m3 <- polr(as.factor(demsatis) ~ east*treat, data = d, Hess=TRUE, method="logistic")
m4 <- polr(as.factor(demsatis) ~ east*treat, data = d, Hess=TRUE, method="probit")


### Stargazer output
stargazer(m1,m2,m3,m4,
          omit.stat = c("rsq", "f", "ser"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(8),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),
          covariate.labels=c("East Germany",
                             "Economic insecurity",
                             "Political underrepresentation",
                             "East Germany * Economic insecurity",
                             "East Germany * Political underrepresentation"), 
          digits = 3)




# Table A9 ----------------------------------------------------------------

### Regressions
m5 <- polr(as.factor(demsatis) ~ treat_applicable, data=d, Hess=TRUE, method="logistic")
m6 <- polr(as.factor(demsatis) ~ treat_applicable, data=d, Hess=TRUE, method="probit")

m7 <- polr(as.factor(demsatis) ~ treat_applicable, data=d, subset=east==0, Hess=TRUE, method="logistic")
m8 <- polr(as.factor(demsatis) ~ treat_applicable, data=d, subset=east==0, Hess=TRUE, method="probit")

m9 <- polr(as.factor(demsatis) ~ treat_applicable, data=d, subset=east==1, Hess=TRUE, method="logistic")
m10 <- polr(as.factor(demsatis) ~ treat_applicable, data=d, subset=east==1, Hess=TRUE, method="probit")


### Stargazer output
stargazer(m5,m6,m7,m8,m9,m10,
          omit.stat = c("rsq", "f", "ser"),
          df = F,
          notes.align="c",
          no.space=T,
          font.size="tiny",
          model.names=F,
          dep.var.labels.include = T,
          dep.var.labels  = "Democracy satisfaction",
          column.separate = c(8),
          column.labels = c("(1=very dissatisfied; 6=very satisfied)"),
          covariate.labels=c("Economic insecurity -- no relevant experience",
                             "Economic insecurity -- relevant experience",
                             "Political underrepresentation -- no relevant experience",
                             "Political underrepresentation -- relevant experience"), 
          digits = 3) 



